function [stt,etamat,smooth_st,yfor,yferr,logL]=kfilterCumSplit(parvec,funcmod,Y,trainvec,solveopt,addsol)
% ====================================================================
% ====================================================================
%% Model solution, second sample stored in structure second (second sample) 
[G,R,C,eu,SDX,Z,structOne,ssvec,structTwo]=feval(funcmod,parvec,solveopt,addsol);

% A. Non-existence in either sample
if isequal(eu,[1;1])==0 || field.euok==0 || isequal(structOne.eu,[1;1])==0
    return
end

% =========================================================================
%% 2. Set-up extended states space 
%     See buildState.m by AB for an alternative 

%% 2.1 Dimensions and time invariant matrices 
[T,ny] = size(Y); 
ns      = size(G,1);
nx   = size(SDX,2);

if size(Z,1)~=structOne.nzHF;error('Rows in Z must match number of variables observed in High frequency'); end 
if size(structOne.Z,1)~=structOne.nzHF;error('Rows in Z2 must match number of variables observed in High frequency'); end 


Z1tilda  = [Z zeros(structOne.nzHF,structOne.nzLF); ...
    zeros(structOne.nzLF,ns) eye(structOne.nzLF)];
Z2tilda =  [structOne.Z zeros(structOne.nzHF,structOne.nzLF);...
    zeros(structOne.nzLF,ns) eye(structOne.nzLF)];

% Constrant vector 
C1   = [C; structOne.A*C];
C2 = [structOne.C; structOne.A*structOne.C];

% casevec: vector determining the particular page of the T and R time 
% varying matrices
casevec =kron( ones(T/structOne.NUMIT,1) , (1:structOne.NUMIT)' );

% type of temporal aggregation
if structOne.flag_tagsum == 1
    sigma =[0 ones(1,structOne.NUMIT-1)];
else
    sigma =(1:structOne.NUMIT);
end


%% 2.2 Page matrices with pages corresponding to quarters {1,2,3,4} 
Tpag   =  zeros(ns+structOne.nzLF,ns+structOne.nzLF,structOne.NUMIT);
Tpag2  =  zeros(ns+structOne.nzLF,ns+structOne.nzLF,structOne.NUMIT);
Rpag   =  zeros(ns+structOne.nzLF,  nx    ,structOne.NUMIT);
Rpag2  =  zeros(ns+structOne.nzLF, nx    ,structOne.NUMIT);
for ii=1:structOne.NUMIT 
    if structOne.flag_tagsum == 1
        Tpag(:,:,ii) = [G zeros(ns,structOne.nzLF);...
            structOne.A*G eye(structOne.nzLF)*sigma(ii)];
        Rpag(:,:,ii) = [R; structOne.A*R];
        Tpag2(:,:,ii) = [structOne.G zeros(ns,structOne.nzLF);...
            structOne.A*structOne.G eye(structOne.nzLF)*sigma(ii)];
        Rpag2(:,:,ii) = [structOne.R; structOne.A*structOne.R];
    else
        Tpag(:,:,ii) = [G zeros(ns,structOne.nzLF); ...
            (1/sigma(ii))*structOne.A*G eye(structOne.nzLF)*(sigma(ii)-1)/sigma(ii)];
        Rpag(:,:,ii) = [R; (1/sigma(ii))*structOne.A*R];
        Tpag2(:,:,ii) = [structOne.G zeros(ns,structOne.nzLF);...
            (1/sigma(ii))*structOne.A*structOne.G ...
            eye(structOne.nzLF)*(sigma(ii)-1)/sigma(ii)];
        Rpag2(:,:,ii) = [structOne.R; ...
            (1/sigma(ii))*structOne.A*structOne.R];
    end 
end 


%% Demeaning using the split sample
%  First Sample (C1)
if any(C1~=0)==1
    Y(1:addsol.rowBegSecond-1,:)=...
        Y(1:addsol.rowBegSecond-1,:)-repmat((Z1tilda*C1)',[addsol.rowBegSecond-1 1]);
end 
% Second sample (second.second.C)
if any(C2~=0)==1;
    Y(addsol.rowBegSecond:end,:)=...
        Y(addsol.rowBegSecond:end,:)-repmat((Z2tilda*C2)',[(T-addsol.rowBegSecond+1) 1]);
end 
Y=Y';

%% Initialization 
Pzero=feval(@lyapunov_symm,Tpag(:,:,casevec(1)),...
    Rpag(:,:,casevec(1))*(SDX')*((Rpag(:,:,casevec(1))*(SDX'))'));

%% 3. Forward filter 

%% 3.1. Dimensions and storage    

vt=zeros(ny,T);
finvt=zeros(ny,ny,T);
kpartg=zeros(ns+structOne.nzLF,ny,T);
logLnc=zeros(T,1); 
yfor =zeros(ny,T);

% Matrices with one additional entry (initialization)
% to recover observables 
stt=zeros(ns+structOne.nzLF,T+1);
ptt=zeros(ns+structOne.nzLF,ns+structOne.nzLF,T+1);
ptt(:,:,1) = Pzero;

% casevec: vector determining the particular page of the T and R time 
% varying matrices
casevec =kron( ones(T/structOne.NUMIT,1) , (1:structOne.NUMIT)' );

W   =eye(ny);
Zdim      =zeros(T,1); 
mat_obspos=zeros(T,ny);


%% 3.3 Start Forward Filter using KF_DK 
% First Part
for ii=1:addsol.rowBegSecond-1;
    
    % Handling of missing observations
    ytt=Y(:,ii);
    
    % Determine W and position of the NAN
    ind =~isnan(ytt);
    rowt =find(~isnan(ytt));
    
    
    ytt=ytt(ind);
    Zdim(ii)=length(ytt);
    Ztt=W((ind==1),:)*Z1tilda;
    mat_obspos(ii,1:Zdim(ii))=rowt;
    dimt=( 1:Zdim(ii) );
      
    yfor(dimt,ii) =Ztt*(Tpag(:,:,casevec(ii))*stt(:,ii));
    [stt(:,ii+1),ptt(:,:,ii+1),logLnc(ii),vvv,fin,kgain]=...
        kf_dk(ytt,Ztt,stt(:,ii),ptt(:,:,ii),Tpag(:,:, casevec(ii) ),...
        Rpag(:,:,casevec(ii))*SDX' );
    
    % Output dimension depends on DIMT 
    vt(dimt,ii)=vvv; 
    finvt(dimt,dimt,ii)=fin; 
    kpartg(:,dimt,ii)=kgain;     
    
    
end
% Run through filter for second part of the sample
for ii=addsol.rowBegSecond:T;
    
    % Handling of missing observations
    ytt=Y(:,ii);
    
    % Determine W and position of the NAN
    ind =~isnan(ytt);
    rowt =find(~isnan(ytt));
    
    
    ytt=ytt(ind);
    Zdim(ii)=length(ytt);
    Ztt=W((ind==1),:)*Z2tilda;
    mat_obspos(ii,1:Zdim(ii))=rowt;
    dimt=( 1:Zdim(ii) );
    
    
    yfor(dimt,ii) = Ztt*stt(:,ii);
    [stt(:,ii+1),ptt(:,:,ii+1),logLnc(ii),vvv,fin,kgain]=...
        kf_dk(ytt,Ztt,stt(:,ii),ptt(:,:,ii),Tpag2(:,:, casevec(ii) ),...
        Rpag2(:,:,casevec(ii))*structOne.SDX' );
    
    % Output dimension depends on DIMT 
    vt(dimt,ii)=vvv; 
    finvt(dimt,dimt,ii)=fin; 
    kpartg(:,dimt,ii)=kgain;
    
    
end

%% 4. Likelihood with Integration Constant  
consInt=-0.5*sum(Zdim)*log(2*pi); 
logLnc=logLnc(trainvec(1):trainvec(2)); 
logL=consInt+sum(logLnc);

%% 5. Truncate filters and obtain initial observations 
yferr=vt';
yfor =yfor'; 
stt=stt(:,2:end); 
ptt=ptt(:,:,2:end); 

%% 6. Disturbance smoother with TV matrices 
% Obtain the Innovations using a disturbance smoother 
etamat=zeros(nx,T); 
smooth_st=zeros(ns+structOne.nzLF,T); 
rmat=zeros(ns+structOne.nzLF,T); 

%% 6.1 Initialize RSTAR & start at t=Nobs 
rstar=zeros(ns+structOne.nzLF,1);
[rstar,etamat(:,end)]=smoothdis(rstar,(structOne.SDX')*structOne.SDX,...
    (Rpag2(:,:,casevec(end))'),(Ztt'),finvt(dimt,dimt,end),zeros(nstilda),vt(dimt,end));
smooth_st(:,ii)=stt(:,end)+ptt(:,:,end)*rstar; 
rmat(:,end)=rstar; 

%% 6.2 Begin Backward recursion 
Ztr=Z2tilda'; 
Gf = Tpag2;
Rf = Rpag2;
Qf=(structOne.SDX')*structOne.SDX;
for ii=(T-1):-1:1;
    % Recall addsol.rowBegSecond is the beginning of the Second Sample 
    % This break must occur exactly when second sample begins
    if ii==addsol.rowBegSecond-1;
        Ztr=Z1tilda';
        Qf=(SDX'*SDX); 
        Rf = Rpag;
    end
    
    % Varying dimension in the backward recurions
    dimt=( 1:Zdim(ii) );
    % [Ztt]'= [Wtt*Z]' = = Z'*Wtt'
    Ztt=Ztr*( W( mat_obspos(ii,1:Zdim(ii)),:)');
    
    %% NOTE on backward filtering 
    %
    % When dealing with time varying G and R matrices need plug 
    % $$ G_{t+1} $$ and $$ R_{t} $$. Note the difference in timing. 
    %
    [rstar,etamat(:,ii)]=smoothdis(rstar,Qf,(Rf(:,:,casevec(ii))'),...
        Ztt,finvt(dimt,dimt,ii),...
        ((Gf(:,:, casevec(ii+1))-Gf(:,:, casevec(ii+1))*kpartg(:,dimt,ii)*Ztt')'),vt(dimt,ii));
    
    smooth_st(:,ii)=stt(:,ii)+ptt(:,:,ii)*rstar;
    rmat(:,ii)=rstar; 
    
    % Recall addsol.rowBegSecond -1 is the end of the First Sample
    % This break must occur exactly when first sample ends
    if ii==addsol.rowBegSecond-1; 
        Gf=Tpag;    
    end 
     
end
a0 = P0*Tpag(:,:,casevec(1))'*rstar;   % Note: this is only correct in the case where a0 = zeros(ns,1) 
etamat=(etamat)';
smooth_st=(smooth_st)'; 
stt  =stt'; 

%% 7. Check Smoother 
% Check that Smooth States are identical if using disturbance smoother (above) vs. state smoother (below) 
% and if can also recover the observables
tol=1e-5; 
% State at time zero give by GG1*azero+Pzero*rstar;
[yCheck,smoothCheck]=kfilterRegSplitSimulation(a0,etamat'); 
disp('Max Discrepancy Smooth vs. Actual Data'); 
maxdifY=comparemat(yCheck,Y'); 
disp('Max Discrepancy State and Innovation Smoother'); 
maxdifS=comparemat(smoothCheck,smooth_st); 
if maxdifY > tol || maxdifS > tol 
    error('Smoother discrepancy exceeds tolerance') 
end 

%% 8. Initial Variance and State per shock (used below for the counterfactual decompositions)
% vInitialPershock: Initial Variance, decomposed per shock 
% sOnePerShock, State at time 1 decomposed, decomposed per shock
vInitialPerShock=zeros(ns,ns,nx);
sZeroSmoothPerShock=zeros(ns,nx); 
sOneSmoothPerShock =zeros(ns,nx); 
for ii=1:nx
    vInitialPerShock(:,:,ii)=lyapunov_symm(Tpag(:,:,casevec(1)),...
        Rpag(:,ii,casevec(1))*SDX(ii,ii)*SDX(ii,ii)*(Rpag(:,ii,casevec(1))'));
    sZeroSmoothPerShock(:,ii) = vInitialPerShock(:,:,ii)*Tpag(:,:,casevec(1))'*rstar;
    etatemp = zeros(nx,1);
    etatemp(ii) = etamat(1,ii);
    sOneSmoothPerShock(:,ii)=Tpag(:,:,casevec(1))*sZeroSmoothPerShock(:,ii)...
        +Rpag(:,:,casevec(1))*etatemp;
end
maxdifSzero=comparemat(sum(sZeroSmoothPerShock,2),a0); 
maxdifSone =comparemat(sum(sOneSmoothPerShock,2),smooth_st(1,:)' ); 
if maxdifSzero > tol; error('Cannot decompose sZeroPerShock'); end 
if maxdifSone > tol; error('Cannot decompose sOnePerShock'); end


%% 9. Counterfactual Decomposition 
% Generate states by feeding each shock at a time, ensure that it recovers
% the original state 
disp('Begin Counterfactual')
countStates=zeros(T,ns,nx);
countObs   =zeros(T,ny,nx);
for ii=1:nx 
    etaTemp=zeros(T,nx); 
    etaTemp(:,ii)=etamat(:,ii); 
    [countObs(:,:,ii),countStates(:,:,ii)]=kfilterRegSplitSimulation(sZeroSmoothPerShock(:,ii),etaTemp');
end 
disp('Maximum Discrepancy Counterfactual States and Smooth States') 
maxdifCount=comparemat(sum(countStates,3),smooth_st); 
if maxdifCount > tol; error('Counterfactual States do not recover smooth states'); end 
   
%% Subroutine kfilterRegSplitSimulation allows to simulate the model  
% Inputs 
% sInitial: [ns 1] Initial State 
% innovMat: [nx T] matrix of innovations 
% By being a sub-routine it has access to all variables defined above 
% Be-careful not to use an index (ii,jj) used above or to repeat variable
% names
    function [ySim,sSim]=kfilterRegSplitSimulation(sInitial,innovMat) 
        if ~isequal(size(innovMat),[nx T]) 
            error('Input innovMat must be [nx T]') 
        end 
        sSim=zeros(ns,T);
        ySim=zeros(size(Y));
        sSim(:,1)=Tpag(:,:,casevec(1))*sInitial + ...
            Rpag(:,:,casevec(1))*innovMat(:,1);
        ySim(:,1)=Z1tilda*sSim(:,1);
        for kk=2:addsol.rowBegSecond-1;
            sSim(:,kk)=Tpag(:,:,casevec(kk))*sSim(:,kk-1)+...
                Rpag(:,:,casevec(kk))*innovMat(:,kk);
            ySim(:,kk)=Z1tilda*sSim(:,kk);
        end
        for kk = addsol.rowBegSecond:T
            sSim(:,kk)=Tpag2(:,:,casevec(kk))*sSim(:,kk-1)+...
                Rpag2(:,:,casevec(kk))*innovMat(:,kk);
            ySim(:,kk)=Z2tilda*sSim(:,kk);
        end
        ySim=ySim'; 
        sSim=sSim'; 
    end
%% End of File
 end 